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^— U ', In this paper we present a dynamical Monte Carlo algorithm which is 

^N I applicable to systems satisfying a clustering condition: during the dynamical 

evolution the system is mostly trapped in deep local minima (as happens in 
glasses, pinning problems etc.). We compare the algorithm to the usual Monte 
Carlo algorithm, using as an example the Bernasconi model. In this model, 
a straightforward implementation of the algorithm gives an improvement of 
several orders of magnitude in computational speed with respect to a recent, 
already very efficient, implementation of the algorithm of Bortz, Kalos and 
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Direct dynamical simulations have traditionally played an important role in statisti- 
cal physics. Investigations have covered a large range of problems, from dynamical phase 
transitions (e.g. the spinodal decomposition and the glass transition) to detailed, realistic 
simulations of polymers and long proteins in the sub-nanosecond regime. 

The macroscopic evolution of the system is always much slower than its microscopic 
Monte Carlo dynamics. This simply means that the problem of dynamical simulations is, 
almost by definition, a very difficult enterprise. 

In dynamical Monte Carlo simulations one considers Markov sequences of configurations 

o"i(t-i) -^ o"2(^2) • • • ^ CTzin) -^ cTi+i{ri+i) . . . ctc (1) 

where the symbols mean that the system will remain in the configuration ai during the 
time Tj. For concreteness we imagine the configuration to stem from a spin model a = 
{Si, 5*2, ... , Sn) with Si = ±1. In the single-spin fiip Metropolis dynamics, the system can 
move from configurations ctj to any configuration which differs from it by the flip of a single 
spin m, which we denote by a^^\ The probability of making a move is then given by 

^^^ ' > iV \exp[-/?(E(aH)-E(a))] [otherwise) ^'••-^^ l^J 

where E[a) is the energy of the conflguration a. Usually, the sequence in eq. (p is directly 
simulated (with a finite timestep Ar) by means of a rejection method. In that case, the 
time intervals r, in eq. (|l]) are multiples of Ar: Tj = M Ar. At each time step a flip of a 
randomly chosen spin m is attempted, and a uniform random number is drawn. Then 

, . . r cr'"^! if pia ^>- a^"^^) > ran < ran < 1 ,„. 

cT r + Ar = { , •^ . ^\ ' (3) 

o [otherwise) 

At low temperatures, the MC dynamics often encounters the problem of small acceptance 
probabilities: as soon as the temperature is such that the available phase space is dominated 
by low-lying minima, the system will take a long time to move from a low-lying state to an 
excited one. The problem of small acceptance probabilities is certainly a nuisance, but not 
a true impediment to the simulation. Almost 20 years ago, Bortz, Kalos and Lebowitz |1| 
(BKL) showed how to circumvent it: in a nutshell, their algorithm consists in giving up the 
rejection method in favor of a direct calculation of waiting times, which can be calculated 
from the knowledge of all the transitions. In the above example, the time r^ the system 
will stay in a conflguration o"j can be directly sampled from its probability distribution 
p(tj) = Aexp(— Arj) with A = Yl,m.=\,NViPi ~^ ^i\ The BKL algorithm has been used 
repeatedly [e.g. in f^ |^) over the last years. A recent, very efficient implementation 0] 
has given improvements in computational speed of the order of 10^ to 10^ for glassy systems 
at low temperatures. 

If the problem of small acceptance probabilities has thus been effectively solved by the 
BKL algorithm, it is usually accompanied by another one, which we may call "futile" dy- 
namics: if the dynamical evolution is dominated by isolated deep local minima, the system 
will get trapped even though we force it to make transitions out of them: very shortly after 
having escaped a minimal configuration (by means of the BKL algorithm), the system sim- 
ply falls back into it. This means that, at low temperature, the dynamics consists mostly of 
short cycles. 



The purpose of the present paper is to show that this problem, also, can be eliminated, 
and the power of the BKL algorithm increased by at least a few orders of magnitude. To 
give an illustration of the problem we are concerned with and to introduce the model which 
will serve as an example for our presentation, we show in fig. 1 the result of a very long 
simulation of a single run of the Bernasconi model |^, defined by the hamiltonian: 
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The model defined in eq. (^) has acquired some fame during the last years as an example 
of a "tough" optimization problem. An inkling of why this is so can immediately be gotten 
from fig. 1, in which we show the temporal evolution of a system of size N = 400 at 
a temperature T = 0.05. The simulation has gone over a time t = 3.9 x 10^, {i.e. has 
performed the equivalent of 3.9 x 10® Monte Carlo steps per spin). The total number of 
accepted moves is around Nacc = 2.1 x 10^, but the total number of different configurations 
visited, N^iff = 1290, is very small. In an effort to be completely explicit, we restate that 
the small acceptance probability {Nacc/{N * t) -C 1) is related but not identical to the 
"futility" expressed by N^iff/Nacc ^ 1- On the second curve in fig. 1 we display the total 
number of different configurations visited as a function of time. One can clearly see that 
this number grows very slowly whenever the system is within a given plateau of the energy. 
Many of the configurations are thus visited a large number of times. It should be evident 
that the total time the system spends in any of these configurations can be evaluated in a 
much more efficient way than by visiting them repeatedly. It can in fact be deduced from 
the Boltzmann distribution. We have extensively tested this fundamental hypothesis of our 
method in simulations using the BKL algorithm. 

Consider now the schematic view of the algorithm presented in fig. 2. This algorithm 
completely avoids the computational slowdown associated with repeated visits to the same 
configuration. In fact, it visits one new configuration per iteration. At any given time, there 
is an active set of configurations which make up the "cluster" C. For any configuration i 
in C, we have stored, besides the energy E(i), the total probability to move from site i to a 
configuration outside the cluster 

Pext{i) = Y.P(^^ 3) (5) 

Now we suppose that the configurations in C are in quasi-equilibrium, i.e. we assume that 
the system visits configuration j with the probability given by its normalized Boltzmann 
weight 

exp(-/?E(^)) 
^^"^^ E.ecexp(-/3E(A:)) ^'^ 

The (normalized) probability to move from any configuration in C to a configuration outside 
C is thus given by Y.ieC Veqij^Pextii) ■ In direct analogy to the BKL algorithm, we then sample 
the time r the system will stay in C from its probability distribution 

p(r) = Aexp(-Ar) with \ = ^Peq{i)Pext{i) (7) 



As before, A is the parameter of the exponential distribution of waiting times in the cluster. 
It should be evident that A can very quickly become a small quantity, which needs to be 
calculated with high-precision arithmetic. In the calculation presented later on (fig. 4), A~^ 
routinely takes on values of A~^ ~ 3 x 10^, which simply means that a single iteration of 
the algorithm (at very large times of the simulation) corresponds to 3 x 10^ Monte Carlo 
steps per spin! After this time, sampled from eq. (|^, the system moves outside C. Now, 
in analogy with the BKL algorithm, we are able to sample in its turn the configuration of 
C (denoted by A in fig. 2) from which this move will be made (this probability is simply 
proportional to Peq{j)Pext{j)) ■ Finally, having identified the configuration ctj (= A), we can 
sample the configuration outside C, according to eq. (0), after having excluded all the moves 
which stay within C. This then produces the configuration denoted by B in fig. 2. 

In our current implementation of the algorithm, we then calculate the probability Pduster 
to return to the cluster from B, or to move directly to configurations outside the cluster, 
Pext- In both cases, the amount of time spent on B is sampled from p^eif = 1 —Pext — Pduster- 
We then add B to the cluster and update all the probabilities Pexti}) for all members of C, 
whose size, in the meantime, has been incremented. 

A little thought shows that the algorithm just sketched can be implemented recursively, 
and that a single iteration consists in two passes through the cluster. Furthermore, besides 
the configurations (which we keep in external memory), we only need the tables E{i) and 
Pext{i), i-G- there is no need for tables which grow faster than the cluster size. In the example 
treated, we can presently, for A^ = 400, allow up to 10^ members in C. 

To show that the algorithm works, we have tested it against the BKL algorithm. In 
fact, for both algorithms, we have started 300 samples on the same configuration (marked 
"X" in fig. 1), and have determined the histogram of trapping times (times after which the 
system accesses a new plateau with lower energy). This histogram is shown in fig. 3. As 
can be clearly seen, the distributions of waiting times are identical, which indicates that the 
algorithm is correct for all intents and purposes. Furthermore, it should be evident that, 
once an old plateau has been left, and a new configuration of minimal energy found, the 
original members of C drop out of the picture. At this point, a new cluster can be started. 

In fig. 4 we show the result of a simulation for the same system as in fig. 1, but this 
time with the new algorithm. Using the same amount of computer time, we reach physical 
times which are about 10000 times larger. The inset in fig. 4 shows the energy of the new 
configurations which enter the cluster. This energy drops as a new plateau is found, and then 
increases (approximately logarithmically in time), as more and more remote configurations 
are encountered during the thermal exploration process. The main plot in fig. 4 shows 
the mean energy of the cluster. Note that many more quantities are accessible, such as 
correlation functions, or more detailed properties of the dynamical process. 

Finally we would like to discuss whether the algorithm, as presented, and as sketched 
in fig. 2 can be made rigorous. In the present state, we increase the size of the cluster by 
one at any iteration. Clearly, a more rigorous approach would be to add point "B" to the 
cluster only if p^xt "^ Pint, say if Pext/Pint < 7- The algorithm is equivalent to the BKL 
algorithm for 7 = 0, but certainly retains some asymptotic validity as 7 > 0. The issue at 
hand is that the configurations in the cluster should satisfy a quasi-equilibrium condition, 
which basically means pext ^ Pint, and which may be implemented with various degrees of 
rigor. 



A second question of course concerns the criterium for abandoning the cluster and grow- 
ing a new one. In our present case we have simply adopted an energy criterium: each 
time a new lowest energy was found, the old cluster was abandoned, and a new one grown. 
The precise criterium will depend on the particular problem. Whenever geometry plays an 
important role, the cluster may have to be scrapped as soon as the geometrical distances 
between members of the cluster become too large. 

In conclusion, we have presented in this paper an algorithm for the Monte Carlo dynamics 
of glassy systems, which is able to produce enormous gains in efficiency if configurations are 
visited over and over again. In a natural way, the algorithm forms a cluster of a certain 
number of configurations, which are close in energy. These configurations are then assumed 
to be in thermal equilibrium, an assumption which was very well satisfied in the example 
we treated. 

Finally, we indicate that a FORTRAN source code for the algorithms used in this work 
may be obtained by email. 
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Figure Captions 

1. Energy vs logt for the Bernasconi model at T = 0.05 for a single realization. The 
calculation was performed using the BKL algorithm in a very efficient implementation 
[^]. The point marked "X" serves as a starting point for subsequent simulations (c/fig. 
3). The second curve gives the total number of different configurations visited (scale 
to the right). When on a "plateau", the number of configurations increases roughly 
linearly in log t. 

2. Schematic view of the Cluster Monte Carlo algorithm presented in this paper. 

3. Histograms of trapping times for the two algorithms. For this figure, both the BKL 
algorithm and the cluster algorithm presented in this paper were started at the point 
marked "X" in fig. 1. 300 runs of both algorithms were done, and the time after which 
a new plateau of lower energy was found was typically between t = 10"^ and t = 10^. 

4. Mean energy of the cluster vs logt for the Bernasconi model at T = 0.05 (single run 
of the cluster algorithm). The inset displays the energy of new configurations which 
enter the cluster. 
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KRAUTH, PLUCHERY : A rapid Monte Carlo. 



